---
title: 'Latta et al 2017 PeerJ: Multilevel Regression to calcualte species-specific slopes using BLUPs'
author: "brouwern@gmail.com"
date: "December 20, 2016"
output:
  word_document: default
  html_document: default
  fig_width: 6
  fig_height: 4
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)
```


Workspce saved as "BLUPS_model_workspace_Jan8_2017.RData"

# Part Ib: Running Multilevel models for BLUPS

* Multilevel regression of Costa Rica bird abundance ~ 1  
* general form: lmer(abund ~ 1 + (1|year) + (time|species))

Goal is to get BLUPS to describe individual-level trends

## Task done in this script

* Load long data
* Run models w/ranefs for year and spp

# Preliminaries

## Load data

```{r, message=FALSE, warning=FALSE}

#setwd
my.dir <-  "C:/Users/lisanjie2/Dropbox/Latta_et_al_2017_PeerJ_FINAL"

setwd(my.dir)

#load data
CR.long <- read.csv("./DATA/Latta_et_al_mistnet_long.csv")

CR.wide <- read.csv("./DATA/Latta_et_al_mistnet_wide.csv")


#load cleaned data in "long" format for main modelling
#load("./cleaned_data/CR_data_long_update12202016.RData")





```


## Load packages

```{r, message=FALSE, warning=FALSE}
library(lme4)
library(arm)  #for getting SEs for BLUPS *
library(ggplot2)

# * see figuraleffect.wordpress.com/2011/07/01/computing-standard-errors-ses-and-confindence-intervals-for-blups-in-lmerlme4/
# for info on SEs of BLUPS

```




# Run hierachical models

## Define indices for subsets

Use only speices w/ < 7 years total not observed (that is, > 3 years observed)

* Jan only
* Aug only
* both months

```{r}

### Define indices for subsets
i.Jan.obs3.or.more <- with(CR.long, 
                           which(month == "Jan" &
                                   yrs.not.obs < 7))
length(i.Jan.obs3.or.more)#680

i.Aug.obs3.or.more <- with(CR.long, 
                           which(month == "Aug" &
                                   yrs.not.obs < 7))
length(i.Aug.obs3.or.more)#950



i.BOTHMONTHS.obs3.or.more <- with(CR.long, 
                           which(yrs.not.obs < 7))
length(i.BOTHMONTHS.obs3.or.more)#950


```



## Models



```{r, cache = TRUE, warning=F}

#model: single model w/both months
glmer.poisnorm.both <- glmer(N ~ month +  1 +
                        (1|name) +               #spp-level intercept
                        (year.cent|name:month) + #slopes for each spp w/in each motnh
                        (1|year) +               #intercept for each year
                        (1|year:month) +         #intercept for each year nested w/in month
                        (1|i),                   #ind level raned for pois-norm
                        family = poisson,
                       data = CR.long[i.BOTHMONTHS.obs3.or.more,])
```




<br>



## Save output of models

### Save model lists

Save output of each model to /model_output folder
```{r}
# save(glmer.poisnorm.jan, file ="./model_output/BLUP_model_Jan_pois_norm.RData")
# 
# save(glmer.poisnorm.aug, file = "./model_output/BLUP_model_Aug_pois_norm.RData")

save(glmer.poisnorm.both, file = "glmer_poisnorm_both.RData")

```



### Function to extract BLUP info from models

Function extact random effects info and does some labeling.

```{r}

### Function do extract information
extract.etc <- function(model,
                        ranef.use = "name",
                        subset){
  
  #betas (BLUPS)
  out.beta <- data.frame(ranef(model)[[ranef.use]])
  
  #SEs
  out.se   <- data.frame(se.coef(model)[[ranef.use]])
  
  #rename
  names(out.beta) <- c("int","slope")
  names(out.se) <- c("int.SE","slope.SE")
  
  #Calc CI
  out.se$slope.95CI <- out.se$slope.SE*1.96
  
  #combine
  out <- cbind(out.beta,out.se)
  
  #make stars
  out$sig <- ifelse(abs(out$slope) > out$slope.95CI,"*","")
  out$sig <- ifelse(abs(out$slope) < out$slope.95CI & 
                      abs(out$slope) > out$slope.SE,".",                   out$sig )
  
  
  #add meta data
  out$subset <- subset
  out$spp <- factor(row.names(out))
  row.names(out) <- NULL
  
  names(out) <- paste("glmer",
                      names(out),
                      subset,
                      sep = ".")
  
  return(out)
}


```

<br>

## Extract BLUP information from ranefs

Use extract.etc() from previous to process data

* Jan only model
* Aug only model
* Both months in same model

```{r}

#Extract "both data
both.fin <- extract.etc(glmer.poisnorm.both,
                       ranef.use = "name:month",
                       subset = "both")

```

<br>

### Do some labeling and clean things up

```{r}

#labels months
both.fin$month <- NA
both.fin$month[grep("Jan",both.fin$glmer.spp.both)] <- "Jan"
both.fin$month[grep("Aug",both.fin$glmer.spp.both)] <- "Aug"




#clean spp names to remove month 
names(both.fin)[8] <- "name"
both.fin$month <- factor(both.fin$month)
both.fin$name <- gsub("[:][AJ][ua][gn]","",both.fin$name)

```


### Subset by month

This subsets each month from the "both months" data.  The results can then be placed side by side.  That is, this essentilly reshapes the data from long to "wide" format.

```{r}

#subset by month
both.fin.Aug <- subset(both.fin, month == "Aug")
both.fin.Jan <- subset(both.fin, month == "Jan")

#focal columns to use from both dfs
colsuse <- c("name",
             "glmer.slope.both",
             "glmer.slope.SE.both", 
             "glmer.slope.95CI.both",
             "glmer.sig.both")

# isolate focal columns
both.fin.Aug2 <- both.fin.Aug[,colsuse]
both.fin.Jan2 <- both.fin.Jan[,colsuse]

# change names so each subset has unique col names
names(both.fin.Aug2) <- gsub("glmer","Aug",names(both.fin.Aug2))
names(both.fin.Aug2) <- gsub(".both","",names(both.fin.Aug2))

names(both.fin.Jan2) <- gsub("glmer","Jan",names(both.fin.Jan2))
names(both.fin.Jan2) <- gsub(".both","",names(both.fin.Jan2))

```

<br>

### Merge the two months

Each species gets its own row.  Months in seperate columns
```{r}
both.fin2 <- merge(both.fin.Aug2,both.fin.Jan2, all  = T)


```



### Indicate if there is sig or marg for either month
```{r}

i.sig <- with(both.fin2, which(Aug.sig != "" | Jan.sig != ""))

both.fin2$one.sig <- NA
both.fin2[i.sig, "one.sig"] <- "one.sig"

```


### Add species info

```{r}
#screen by yrs observe
i.use <- which(CR.wide$yrs.not.obs < 7)

spp.info <- CR.wide[i.use,c("family",
                    "hab.prime",
                    "forage.guilde",
                    "name")]
i.dup <- which(duplicated(spp.info$name) ==TRUE)
spp.info <- spp.info[-i.dup,]


both.fin3 <- merge(spp.info,
              both.fin2, 
              by = c("name"), all = T)
dim(both.fin2)
dim(both.fin3)


```



### Function to round numbers in a dataframe    
```{r}
# from http://stackoverflow.com/questions/9063889/how-to-round-a-data-frame-in-r-that-contains-some-character-variables
round_df <- function(df, digits) {
  nums <- vapply(df, is.numeric, FUN.VALUE = logical(1))

  df[,nums] <- round(df[,nums], digits = digits)

  (df)
}
```





### Exponentiated
```{r}

exp.slope <- apply(both.fin3[,c("Aug.slope","Jan.slope")], 2, exp)

CI.up.Jan <- with(both.fin3, exp(Jan.slope + Jan.slope.95CI))

CI.lo.Jan <- with(both.fin3, exp(Jan.slope - Jan.slope.95CI))


CI.up.Aug <- with(both.fin3, exp(Aug.slope + Aug.slope.95CI))


CI.lo.Aug <- with(both.fin3, exp(Aug.slope - Aug.slope.95CI))

#paste Jan CIs together
Jan.CIs <- paste("(",round(CI.lo.Jan,2),
      "-",
      round(CI.up.Jan,2),
      ")",
      sep = "")

Jan.CIs <- gsub("\\(NA-NA\\)","",Jan.CIs)


#paste Aug CIs together
Aug.CIs <- paste("(",round(CI.lo.Aug,2),
      "-",
      round(CI.up.Aug,2),
      ")",
      sep = "")

Aug.CIs <- gsub("\\(NA-NA\\)","",Aug.CIs)

Species <- gsub("(.*)([A-Z])","\\1 \\2",
              both.fin3[,c("name")])
Species <- gsub("^ ","",Species)
Species <- gsub("White-breastedWood- Wren",
                "White-breasted Wood-Wren",Species)


Species <- gsub("Orange-billedNightingale- Thrush",
                "Orange-billed Nightingale-Thrush",
                Species)
Species <- gsub("Red-crownedAnt- Tanager",
                "Red-crowned Ant-Tanager",
                Species)




Species2 <- paste(Species, " (",
                  both.fin3$hab.prime, "-",
                  both.fin3$forage.guilde,
                  ")",
                  sep = "")

#Augst slope - round and remove NAs
Aug <- round(exp.slope[,"Aug.slope"],2)
Aug[is.na(Aug)] <- ""

#Jan slope - round and remove NAs
Jan <- round(exp.slope[,"Jan.slope"],2)
Jan[is.na(Jan)] <- ""

#Add sig star
Aug.sig <- both.fin3$Aug.sig
Aug.sig[is.na(Aug.sig)] <-  ""
Aug <- paste(Aug, Aug.sig,sep = "")

Jan.sig <- both.fin3$Jan.sig
Jan.sig[is.na(Jan.sig)] <-  ""
Jan <- paste(Jan, Jan.sig,sep = "")

both.fin4 <- data.frame(Species = Species2,
                   August = Aug,
                   Aug.CI  = Aug.CIs,
                   January = Jan,
                   Jan.CI  = Jan.CIs,
                   both.fin3[,c("one.sig",
                             "family")]
                   )



```


Save all species
```{r}
#write.csv(both.fin4, file = "blups_all_spp.csv")
```




```{r}
#remove nonsig
both.fin5 <- subset(both.fin4, one.sig == "one.sig")

both.fin5 <- both.fin5[order(both.fin5$family), c(1:5)]

```


### save
```{r}

#write.csv(both.fin5, file = "./model_output/BLUPs_output.csv",row.names = F)

```




### Round table
```{r}
both.fin3.round <- round_df(both.fin3, digits = 3)
```

### Display table
```{r}
cols.j <- c("name",
                     "Aug.slope","Aug.sig",
                     "Jan.slope","Jan.sig")
sig.i <- with(both.fin3.round, which(one.sig == "one.sig"))


table.fin <- both.fin3.round[sig.i,cols.j]

table.fin[is.na(table.fin)] <- ""

dim(table.fin)
```



Save table
```{r}
#write.csv(table.fin, file = "blup_table_temp.csv")
```




# Plot Aug. vs. Jan

Guild names
```{r}
both.fin3$forage.guilde <- gsub("F","Frugivore",both.fin3$forage.guilde)

both.fin3$forage.guilde <- gsub("^I","Insectivore",both.fin3$forage.guilde)
both.fin3$forage.guilde <- gsub("^N","Nectarivore",both.fin3$forage.guilde)
both.fin3$forage.guilde <- gsub("^O","Omnivore",both.fin3$forage.guilde)

```

Primary habitat names
```{r}
both.fin3$Habitat <- ifelse(both.fin3$hab.prime == "F","Prim. For","Secondary")
```



```{r, cache=TRUE}
both.fin3$Significant <- ifelse(both.fin3$Aug.sig == "." |
                               both.fin3$Jan.sig == ".", "Marg","0")

both.fin3$Significant <- ifelse(both.fin3$Aug.sig == "*" |
                               both.fin3$Jan.sig == "*", "Sig",both.fin3$Significant)


both.fin3$forage.guilde <- factor(both.fin3$forage.guilde)



both.fin3$name2 <- ifelse(both.fin3$Significant != 0,both.fin3$name,"")


ns <- c("BANQUI","BBGR","BUFG",
"CHETAN","CCTH","GREHER",            
"OBFC","RCAT","RTHN",    
"SBHB","STTA","SLST",              
"SBHB","STHR","SRFC",    
"VASE","WBWW","WRMA",         
"WTTR") 

both.fin3$name2[which(both.fin3$name2 != "")] <- ns

 


```



### Plot Jan vs. Aug comparison on real scale (exp-ed)

Plot on real scale (response scale).
exp(beta)
Model is pure random slopes w/ no fixed effect time trend.
Therefore these are the slopes of the lines for each species

```{r, message=FALSE, warning=FALSE}
plot.out <- ggplot(aes(y = exp(Jan.slope), 
      x = exp(Aug.slope), 
      shape =Significant,
      color = Habitat,
      label = name),
      data = both.fin3) +
  #geom_point(aes(size = Significant)) +
  facet_wrap(~forage.guilde, scales = "free") +
  geom_text(aes(label=name2),size = 6) +
  xlab("August Trend") +
  ylab("January Trend") +
  geom_hline(yintercept = 1) +
  geom_vline(xintercept = 1) +
  geom_errorbar(aes(x =    exp(Aug.slope),
                    ymin = exp(Jan.slope-Jan.slope.SE),
                    ymax = exp(Jan.slope+Jan.slope.SE))) +
  geom_errorbarh(aes(y = exp(Jan.slope),
                    xmin = exp(Aug.slope-Aug.slope.SE),
                    xmax = exp(Aug.slope+Aug.slope.SE))) +
  theme_bw() +
annotate("text", x = 1.1572, y = 1.3, label = "Inc. both months" ) +
annotate("text", x = 0.885, y = 0.8, label = "Dec. both months" ) +
annotate("text", x = 1.1572, y = 0.8, label = "Inc. Aug/Dec. Jan" ) +
annotate("text", x = 0.885, y = 1.3, label = "Inc. Jan/Dec. Aug." ) 
```




### Make coefficients plt for all species

Plot slope for each species
```{r, message=FALSE, warning=FALSE}
#order by coefficient
i.order <- order(both.fin$glmer.slope.both)

both.fin$i[i.order]<- 1:dim(both.fin)[1]
both.fin$i <- both.fin$i



qplot(y = exp(glmer.slope.both),
      x = i,
      data = both.fin
      #shape = lme.sig.both,
      #color = hab.breadth
      ) + 
  geom_errorbar(aes(ymax = exp(glmer.slope.both+glmer.slope.95CI.both),
                    ymin = exp(glmer.slope.both-glmer.slope.95CI.both)))+
  geom_hline(yintercept = exp(0)) + theme_bw()


```


## Plot trends for individual specis



### Make spp-level predictions from model

Keep random slopes and random intecepts for each spp
```{r}

pred.out <- predict(glmer.poisnorm.both, 
         re.form = ~ (1|name) +                #sp int
                              (year.cent|name:month) #+  #sp slp by mo
                              #(1|year) +                #yr int
                              #(1|year:month) +          #yr/mo int 
                              #(1|i)                     #for pois-norm
        )
        
```

Make df of predictions and orig data
```{r}
pred.out2 <- data.frame(glmer.poisnorm.both@frame[,c("N","month","name","year")],
                   pred.out)


head(pred.out2)
```


Put on real scale
```{r}
pred.out2$Predictions <- exp(pred.out2$pred.out)
```

Focal species from table.fin
```{r}
i.use <- which(pred.out2$name %in% table.fin$name)

```

Add species codes for focal species
```{r}
spp.codes <- read.csv("./DATA/Latta_et_al_species_codes_for_spp_with_sig_trends.csv")

spp.slopes.df <- pred.out2[i.use,]
spp.slopes.df$group <- with(spp.slopes.df, paste(name,month))


spp.slopes.df <- merge(spp.slopes.df, spp.codes)


spp.codes$code



spp.slopes.df$code2 <- factor(spp.slopes.df$code,
       levels = spp.codes$code[order(spp.codes$order)])

```



### Plot slopes for each species
```{r}
names(spp.slopes.df)[3] <- "Month"

base.out <- ggplot(data = spp.slopes.df,
       aes(y = Predictions,
           x = year, 
           color = Month,
           group = group)) +
  geom_line(aes(linetype = Month)) +
  geom_point(aes(y = N, x= year)) +
  facet_wrap(~code2,scales = "free") +
  theme_bw() +
  theme(strip.text = element_text(size = 8, 
                                    face = "bold")) +
  theme(strip.background = element_rect(fill="white")) +
  xlab("Year") +
  ylab("Net Captures per year")



```


Finalize
```{r}

plot.out <- base.out +  theme(axis.text=element_text(size=7)) +
   scale_x_continuous(breaks = c(2005,2008, 2011,2013))

  
ggsave(file = "Appendix_S6.png",plot = plot.out,width = 9.691,height = 5.687)
 
```

